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' Abstract 
Oh' 

' We discuss the hyperfine shifts of the Positronium levels in a relativistic framework, starting 

, from a two fermion wave equation where, in addition to the Coulomb potential, the magnetic 

' interaction between spins is described by a Breit term. We write the system of four first order 

I differential equations describing this model. We discuss its mathematical features, mainly 

in relation to possible singularities that may appear at finite values of the radial coordinate. 
^S) ' We solve the boundary value problems both in the singular and non singular cases and we 

^ develop a perturbation scheme, well suited for numerical computations, that allows to calculate 

the hyperfine shifts for any level, according to well established physical arguments that the 
Breit term must be treated at the first perturbative order. We discuss our results, comparing 
, them with the corresponding values obtained from semi- classical expansions. pacs03.65.Pm, 

^ ■ 03.65.Ge 
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! 1 Introduction 



The relativistic description of the fine structure of the hydrogen atom levels was first proposed 
by Darwin T in a semi-classical treatment of the Dirac equation. This was immediately followed 
' by the Breit proposal of a two body relativistic equation 2 that, in addition to the Coulomb 

potential, included a quasi-static magnetic term where the velocities were substituted by the Dirac 
7?-matrices according to a proposal of Heisenberg. Shortly later Fermi calculated the spectrum 
of a Dirac electron interacting with a Pauli nucleus and deduced the values of nuclear magnetic 
momenta from the measured hyperfine splitting 3 , using the Schrodinger non relativistic wave 
functions to calculate the averages near the origin. In fact, following a Pauli suggestion, also 
the Breit paper contained a semi-classical expansion along the lines of the Darwin work, but 
eventually Breit himself was doubtful about the actual predictivity of his equation and two years 
later he reached the conclusion that the Coulomb interaction should be treated exactly, while the 
magnetic term had to be considered as a first order self-consistent perturbation, so that only its 
diagonal matrix elements on the exact Coulomb eigenstates were brought to bear. Since then, 
major achievements were obtained by Pirenne Berestetski and Landau 5 . Still using the Breit 
semi-classical approximation and the non relativistic Coulomb wave functions, they were able to 
produce analytical approximations for the shifts of the Parapositronium levels and for the ground 
state of the Orthopositronium. In their approach they also included the annihilation term. In the 
first fifties, the appearance of the Schwinger [Sj and of the Bethe-Salpeter equations [71|Sj made 
it possible to calculate up to the order {a = fine structure constant) the shifts of the first 
excited levels of the Positronium, by accounting for the effects of one and two virtual photons, self 
energy and vacuum polarization, but still keeping the semi-classical perturbative scheme [51 ll(J|. 
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Later, in the seventies, due to the great improvement in the experimental analysis of the atomic 
spectra, , to the qualitative changes in the mathematical and physical framework of symmetries 
and mainly to the new ideas of hadrons as composed by quarks, the interest for the bound states 
and for the relativistic wave equations raised up again and never weakened since, |12 | -|24 ) . The 
challenge of the completely relativistic calculation of the hyperfine splitting was again pushed on 
foreground and various models for two spin i interacting particles were proposed, with special 
attention to the Positronium that constituted an ideal system both from a theoretical and an 
experimental point of view: it appears, however, that the use of the semi-classical expansion and 
non relativistic Coulomb wave functions as a starting point has maintained his role. Even in a 
later paper |"131 , where a non perturbative treatment is claimed, the solutions of the Breit equation 
are calculated analytically but expanding the equation up to the second order in a. This latter 
procedure, in particular, may be assumed somewhat safely in atomic physics, where the velocities 
are of the same order of the fine structure constant, so that the expansion in a in fact corresponds 
to a semi-classical approximation. It becomes less and less justified in view of an extension of the 
method to quark bound states (see e.g. [251 126| ). where an expansion in the coupling constant is 
not allowed and where it would be worth dealing from the very beginning with relativistic states. 

In the present paper we fill this old gap and we present a completely relativistic treatment of the 
hyperfine splitting based on the Breit approach, providing an effective method for its computation 
for any spectral level. We thus analyze an approximate interaction within the framework of an 
exact kinematics. It is rather evident that analytic results will be possible only for the initial steps 
of our treatment, i.e. for establishing the system of equations to be discussed and their reduction 
due to conserved quantities. Analytic expressions will also be available when looking for series or 
asymptotic solutions, but finding the spectrum will necessarily be achieved by numerical methods. 
Moreover we shall omit from our treatment the annihilation term. 

The starting point is thus the two fermion relativistic wave equation we presented in [27j . In 
that paper, using a canonical reduction of the relativistic kinematics of the two body problem, we 
introduced Lorentz invariant interactions dependent upon the reduced coordinates and we gave a 
solution of the relative time problem. We then quantized the model assuming a fermion nature 
for both of the two particles and we deduced a completely Lorentz covariant internal dynamics, 
that was reduced by separating the radial part in a multipole scheme exploiting the conservations 
of angular momentum and parity. A further reduction of the radial equations produced a final 
linear system of order four containing the spectral parameter to be determined. We proved that 
the Dirac and the non-rclativistic limits were recovered and we compared our model with other 
existing models, JJ, JX -LSj -12, j finding a general agreement. The complete spectral curves from 
Dirac to Positronium for pure Coulomb interaction were plotted for ground and higher states and 
the crossing of terms inferred in '17' was precised and made concrete. In that paper we also looked 
at the possibility of using the eigenfunctions of the degenerate singlet-triplet ground states of the 
Positronium to set up a direct perturbative calculation of the hyperfine splitting, using the usual 
perturbative terms for the magnetic interaction: the answer was negative even for the ground 
state, since the behavior of the relativistic wave functions in the origin was not compatible with 
the magnetic interaction terms obtained by the semiclassical approximation. 

In order to investigate the hyperfine structure, we shall therefore add to our pure Coulomb 
wave equation a Breit term in a Lorentz covariant way. The radial reduction of the problem is 
again done by using the still holding conservations of angular momentum and parity and here also 
the final wave equation reduces to a fourth order linear system. Unfortunately the situation is 
now less nice than in the pure Coulomb case, since a singularity of the wave equation appears 
at finite positive values of the radial coordinate. This occurrence had already been notified (see 
|28j ) and, to our knowledge, has prevented up to date any correct integration attempt. However a 
careful investigation shows that the behavior of the system is not so bad. Indeed the singularity 
can be "bridged" and the spectrum can be determined although, as previously said, only the 
first perturbative order in the Breit interaction makes sense and higher order corrections should 
be dealt with using the QED, Therefore we are going to present a perturbative approach 

that allows the complete numerical calculation of the contribution to the shift of any level of the 
Para and the Orthopositronium due to the Breit term. We shall show that in almost any case the 
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presence of the additional singularity can be avoided: only for the levels with odd parity and angular 
momentum j = this is impossible and the analysis has to be refined. However, in order to provide 
numerical evidence of the correctness of the Breit argument concerning the perturbative nature of 
the magnetic interaction, we further investigate the singular cases connected to the ground states. 
A numerical approach based on a different starting point and on Pade approximation techniques 
has been successful in calculating the shifts of the Parapositronium singlets to the order a^, |29| . 
The principal merit of our scheme from a numerical point of view is that we are able to obtain the 
results uniquely by means of the numerical values of spectral levels, without any need of using the 
eigenf unctions in order to calculate the matrix elements for the Breit term: in the full relativistic 
case this leads to a faster and more precise approach for finding the spectrum, especially for 
the Orthopositronium triplets. When comparing our results to the values of the hyperfine levels 
computed in the semi-classical scheme and expressed in terms of simple fractions of powers of 
, |1(J[ l3Uj . we find an excellent agreement up to the order a^. This is rather remarkable, since 
the pure Coulomb contribution - where the correct kinematics and the recoil effects are exactly 
accounted for - and the magnetic perturbative terms are here separately different with respect to 
the corresponding semi-classical contributions (more details are given in Section IV). We want to 
stress that, apart some technical obstacles, our general method has a great conceptual simplicity: 
besides its applications to any excited level of electromagnetic bound states with components of 
arbitrary mass ratio, we believe that it should also be relevant for models of different nature, for 
which non-rclativistic calculations are less reliable. 

The paper is organized as follows. In Section II we briefly recall the derivation of the wave 
equation for the two fermion system with a general radial interaction and we present the new equa- 
tion with the additional Breit term together with its reduction that defines the spectral problem to 
be solved. The procedure is completely parallel to that of the pure Coulomb interaction explained 
in the previous paper |27) . which we refer to for detailed expressions, as, for instance, the explicit 
form of the even and odd state vectors. The nature of the mathematical problem and the numerical 
methods we have used for its solution are explained in Section III. Here we also investigate the 
possible presence of singularities at a non vanishing finite value of the radial coordinate and we 
comment on their properties. Finally, in Section IV, we present the results and we discuss them. 
A revisitation of the perturbative expansion as we have used it in our computations (up to second 
order in the proof, but evidently valid at any order) is given in the Appendix. 

In the following we use units such that h = c ~ 1. 



2 The system for the Coulomb and Breit interaction 

In the paper |27| we have described the method for obtaining a Lorentz covariant wave equations 
for two fermions interacting by means of a scalar potential. Obviously we cannot reproduce the 
derivation here and we have to refer to | 27| for details. It is however necessary to recall the main 
definitions in order to have a minimum of self consistency of the treatment. 

(a) The canonical variables. 

Denote by a;^^^ and p|'^^ the Minkowski coordinates and the momenta of two pointlike fermions 
with masses rrii, i — 1 , 2. Define the tensor 



(2.1) 



where = p^^^ + p^^-^ is the total momentum. It satisfies the identities 

^^.e^{P)e1,{P) = ry„/3 e^(P) £;^(P) - ry^" (2.2) 

and therefore it represents a Lorentz transformation to the P = reference frame. We can use 
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(|2.1fl for the construction of a canonical transformation to the variables 

= + ^= ^ + [qar-raq] + qr 

g = f^O 9m ' ^ = ^O'^M ' Qa^^aQt^' ^al^f, (2.3) 



where 



2 (^(1) + ^(2)) ' ^'^—^(1) ^(2)' ^ 2 (^(^) ^(2))' 



Both Ta and Qa are Wigner vectors of spin one, as well cis Z(i is R Newton- Wigner position vector 
for a particle with angular momentum La = Sabc ft Qc ■ In terms of H2.3|l , the two particles mass 
shell conditions p^^^ = mf can be put into the form 

i^qaqa + mlj +yqaqa + mlj =A, X q = -{ml - ml) , (2.5) 

where A — V while the variable q can be fixed and generates a canonical reduction of the phase 
space. Its conjugate coordinate, namely the relative time coordinate r, is cyclic and becomes a 
kind of a gauge function that is chosen a posteriori in order to recover the complete Minkowski 
description for each of the two particles. For instance, a useful choice could be r = 0, although there 
is no necessity of requiring such condition. It is now straightforward introducing a Lorentz covariant 
interaction by changing the mass shell condition H2.5|l with the addition of a scalar potential 
depending upon the Lorentz invariant relative separation r = (rar^)^/^. Thus a relativistic 
two-body system interacting by means of a potential V{r) is described by H2.5|l where A will be 
substituted by h{r) = A — V{r). 

(ft) The two fermion quantum system. 

Let us now quantize the equation (|2.5() assuming that both particles are fermions. In order to 
determine the structure of the wave equation, it is sufficient to consider the case of two free fermions, 
since the interaction will then be introduced in the way previoulsy described. The two particle wave 
function is simply the tensor product of the two single particle wave functions, and must therefore 
satisfy the two separate Dirac equations determined by the operators Di — (^P^-l-fZ^) 7^]^-) —mi and 
D2 — i^Pf^ — qfi) 7(^2) ~"^2, where we adopt the notation 7^^^ — 7''(g)l4, 7^2) = l4®7'^ and where I4 is 
the unity matrix in four dimensions. Introducing the new set of 7-matrices 7 = 7(P) — eQ(P)7^, 
7(j = ^a{P) — ^alfj. using the canonical operators corresponding to the variables (I2.3|) . we 
easily find 

^ = ga(7(i)7(i)<j- 7(2) 7(2) <j) +7(1)^1+7(2)^2, A g = i(m^ - 77?!) . (2.6) 

The system (12.6(1 has exactly the same content as the initial system of the two independent Dirac 
equations. However, put in this form, we see that the variable q remains fixed, in complete 
agreement with the classical canonical reduction. The cyclic character of r is also evident from 
the Lorentz scalar identity for the phase of plane waves P(i)^ ^ -f P(2)'^ a^(2) ^ — P^Z^ — qa^a ■ 
Moreover the definition of 7,7a is actually a unitary transformation on the 7'' 4-vector. Hence, 
as long as P is conserved, the matrices can be represented by the usual 7 matrices. The first of 
equations (|2.6|) is thus the Lorentz-invariant equation for the two-fermion free system. Its sixteen 
eigenvalues are immediately calculated, yielding the expected four singular values 

A = ±{qa qa + ™?)l/2 ± {qa qa + m2)l/2 , ± {q^ + ^2)1/2 ^ ^ ^2)1/2 ^ (2.7) 

each singular value having multiplicity four. 

The equation for the interacting system is simply obtained by substituting A with h{r) = 
A — V{r) in the first of the 12.6|l equations. For future discussions, we find it useful to introduce 
the mass parameters M = mi + m2 , /i — mi — m2 , p — p/M and to reorder the basis of the 
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states making a linear transformation on the tensor product of the two spinor spaces such that 
the system at rest is diagonahzed and the four singular values (|2.7f) are put in the order M, —M, 
—fi and fi. We then make a further linear transformation such that, in each four dimensional 
eigenspace corresponding to the above energy eigenvalues, the square and the third component of 
the total spin S = l4®(T + a^l4 (where a is the Dirac spin) are also diagonahzed and ordered 
with the triplet always following the singlet. 

Wc finally remind that the global parity transformation is given by the product of orbital and 
internal parity transformations. In our picture the internal parity is 7 (gi 7 = diag(l8, —Is)- It can 
be verified that the global angular momentum J = L + S and the parity are conserved, so that, 
together with A, they provide a classification of the states of the global symmetry. 

(c) The radial equations 

In the basis we have chosen, the free Hamiltonian operator Hq is a 16 x 16 matrix of the form 




(2.8) 



where J^a = Adiag(l4, — 14), A = M, /z and T-Ca is a 8 x 8 matrix whose elements are spherical 
differential operators (see for the explicit form). We next construct the "even" and "odd" 
states 'i>± with assigned angular momentum {j,m) and given parity (— )■' or (— )-'^^, whose 16 
components are collected in 4 groups indexed by the eigenvalues ±A of the free system at rest. 



(M) 



(-M) ^(-M) 



(m) 



(2.9) 



In each group the components are singlet-triplet ordered, namely 



(A)^t/^(A) 



± 



(2.10) 



where the subscript "0" refers to the singlet component, while "1+, lo, l-"denote the triplet com- 
ponents. The explicit expressions have been determined using the standard algorithms of the 
composition of angular momenta by means of Clebsch-Gordon coefficients and are therefore linear 
combinations of spherical harmonics. They are fully reported in |27j . 

By applying the Hamiltonian operator (|2.8|l to the states (|2.9|l we obtain a system of radial dif- 
ferential equations by requiring the vanishing of the coefficients of the different spherical harmonics 
in each component of the resulting vector. Indeed we get a large number of differential equations 
{e.g. 34 when starting with ^ — ^'_|_ ), but of course, as one should expect, only eight of them are 
independent for each state with definite parity. Moreover a closer look to them shows that four 
of these eight equations are algebraic relations that can be used, by an appropriate choice of the 
unknown functions, to obtain a system of only four differential equation for each of the two state 
vectors and \E'_. Finally, adding a Lorentz invariant interaction as specified in item (a), the 
fourth order system reads 



dY{r) 
dr 



+ MY{r) =0, 



(2.11) 



where Y{r) — \yi{r), 2/2(^)1 Usi'''), UA^if)) ^i^d is a matrix with general structure 



M 






E{r) 


F{r) 





E{r) 


1 

r 





-F{r) 


G{r) 





2 

r 


E{r) 





-G(r) 


E{r) 


1 



(2.12) 
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The explicit expressions of E{r), F{r) and G(r) for the even case are the foUowing ones 



E{r) 



r h{r) 



F(r) 



2h{r) 



G(r) 



h{r) / + 4 j (j + 1) 



(2.13) 



and they speciahze to the Coulomb interaction when h{r) = A + a/r , a being the fine structure 
constant. 

As explained in |27| . the odd coefficients are obtained from the previous ones by a parity 
transformation, whose action simply results in the change M — /x and n —M. In j^T] it has 
also been shown how the Dirac and the non relativistic limits are recovered from ()2.12|) . We finally 
observe that for E(r) — {e.g. for j = or for /i = in the even case) the fourth order system 
splits into separate second order subsystems, making it easier the numerical solution of the spectral 
problem. In the general case, however, the complete fourth order system has to be considered. 

(d) The addition of the Breit term 

The spin-spin interaction can be described by introducing a Breit term, in addition to the 
Coulomb interaction, in the relation (12.611 that becomes 



A+ - 

r 



f / 

2 V'^(l)^(l)aT'(2)7(2)<j + (7(1)7(1)^ y) (7(2)7(2)b — ) 

= '?a(7(i)7(i),- 7(2)7(2),) +7(i)mi+7(2) 



7712. 



(2.14) 



The Hamiltonian matrix must be modified with respect to H2.8|l in order to account for the presence 
of the additional 7-matrices. This modification as well as the change of the basis are straightfor- 
ward. The angular momentum and parity are conserved and the radial equations are again deduced 
by applying the new hamiltonian to the states 5'±, as in the pure Coulomb case. The general fea- 
tures of the systems of radial equations are also preserved and the final result is again a fourth 
order linear system. In the even case and with h{r) = \ + a/r, the matrix B of this system reads 



B = 





Ee{r) 

GiAr) 




E{r) 

1 

r 



G2,s(r) 



F{r) 



2 

r 

E{r) 





Fe{r) 

EAr) 
1 



(2.15) 



where E(r) and F(r) are given in H2.13|l . The remaining matrix elements read 

(/l2(r)-^2)^2 



Ee{r) = 

GiAr) - 
G2,e(r) = 



r h{r) 
h{r) 



2ae ' 
Aae 



Fe{r) 
4j(.? + l) 



{2as) 



2r{rh{r) - 2ae) 
2ae — rh{r) Aas — rh{r) 



2r 



2j{j + l) Aa^ + {~h{rY + NP)r^ 
h{r) 



2r{rh{r) - 2ae) 



(2.16) 

/i and /i — > — Af. It is also 



Again the matrix of the odd system is obtained by changing M - 
immediate to ver ify that for e = (|TT5|l reduces to fll^ . 

Some remarks are in order. In the first place, we observe that in the matrix (|2.15(l we have 
introduced a parameter e, not present in (|2.14|) and reproducing the latter for the £ = 1/2. We 
shall see in the following that an appropriate use of this parameter permits the calculation of the 
first perturbative terms in the Breit interaction in a way that is numerically much more efficient 
than the usual computations by means of the eigenfunctions. Secondly in H2.14|l and therefore in 
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the system produced by H2.15|) no anomalous magnetic moment is present. This absence would 
constitute a drawback for the calculation of the hyperfine structure of the hydrogen atom, but 
remains an acceptable approximation for the positronium, which is the only case we are going to 
consider later on. Finally, in the even and odd systems describing the positronium, we will let 
/i = and M — 2too, where mo is the electron mass. 



3 The numerical treatment of the Breit interaction 

We now discuss some general properties of the equations for the positronium with a particular 
attention to the possible presence of singularities other than the origin and the infinity. We find 
it useful to introduce the new dimensionless independent variable x and the new eigenvalue w as 
follows 



w 



{^^a^y\\^2m^). (3.1) 



The value of the fine structure constant has been assumed as a = 0.007297372568, ^J. In the 
following we distinguish the discussion of the even case, that develops essentially according to the 
classical O.D.E. theory |S2|, from the odd case, that poses some new problems. 

(a) The even case 

We begin from the case with even parity. From H2.15|l we see that B12 = B21 = B34 = B43 = 
when /X = 0. This means that the fourth order system generated by H2.15|l decouples into two 
second order systems for the unknown functions {yi{x),y3{x)) and (3/2 (a^), 2/4(2;)) that, in turn, can 
be reduced to two second order differential equations for yi{x) and 2/2 (a^)- It has also been shown 
in 1231 that for j = only the system for {yi{x), y^ix)) makes sense, while for j > both systems 
contribute to determining the levels. The first of these two equations, introducing the unknown 
fimction u{x) defined by yi{x) — [(4 + w) x + 2o}^l'^ x"^/^ u(x)^ reads 

u{x) - 



dx"^ 

1 ( 



(4 + w) a; - 2a (4£ - 1) ((4 + w) x + 2a)2 x"^ 
( ~c? ui (8 + w) a;^ — 4 a (4 + w) (2 e + 1) a; - 4a2(4 e + 1) + 16 j (j + 1) ) 
4aej(j + 1) 



((4 + a2 w) a; - 2a {2e- \))x'^ 



u{x) ^ (3.2) 



It can be seen that, in addition to the usual singularities in the origin and at infinity, the equation 
presents a new singularity in the point x = 2a {Ae — l)/(a'^w + 4), that assumes finite positive 
values for e > 1/4. Although we shall show that solutions exist also for e > 1/4 - in fact we shall 
give an explicit solution for e = 1/2 -, in carrying out our perturbative program we can avoid this 
singularity, as well as almost all those we shall encounter in later developments. Indeed, according 
to the perturbative approach we are going to explain here below - and whose proof is given in 
Appendix A -, it is sufficient to solve the equation for values e < 1/4. There is, however, one 
instance in the case with odd parity where the singularity must be explicitly taken into account. 

If 2/2(2;) = [{A + w) X + 2 a {2 e + l)]^^"^ x^^^^ z{x), the second equation of the even case in 
z{x) is 



dx 

1 / 
16a;2 



8ae 3a2(2£+l)2 



{A + a^w)x + 2a{l-2£) {{A + w) x + 2 a {2 e + l))"^ x'^ 
[-a^ w (8 + w) - 4 a (4 + a^ w) (2 £ + 1) a; - 4 (2 e + 1)2 + 16 j (j + 1)) 
4aej(j + l) 2a(2£ + l) 



{{A + w) X + 2 a) x"^ ((4 + w) a; + 2 a (2 e + 1)) a;^ 



z(a:) = 0. (3.3) 



This equation develops a new singularity at a; = 2Q;(2e— 1)/(4 + q;'^w), that assumes positive values 
for e > 1/2. This new singularity is therefore not effective for the same reasons explained above. 
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It would be such also when considering the Breit term as non perturbative, and it would only act 
through a modification of the singularity in the origin. 

The boundary value problem posed by H3.2|l and H3.3|) for sufficiently small values of e, namely 
without the additional singularity, is very classical. Both the boundary points, the origin and the 
infinity, are singular and what we have to do for starting the numerical procedure, is to determine 
the initial conditions in the neighborhood of those points by looking for analytic approximations 
of the solutions in the form of series or asymptotic expansions respectively. We then apply the 
double shooting method in order to determine the value of the spectral parameter with the desired 
accuracy. As a matter of fact both the series and the asymptotic expansions produce two linearly 
independent solutions, but in each case only one solution can be accepted: we are thus in the 
so called 'Himit poinf^ case of the Weyl classification of the boundary value problems [32] and 
therefore there is no ambiguity in choosing the appropriate solution for giving the initial conditions 
and starting the numerical integration, once the spectral parameter has been assigned a numerical 
value that will be adjusted in the successive integrations until when the two solutions coming from 
zero and infinity, as well as their derivatives, match within the required accuracy. This is actually 
the spectral condition that simply reduces to checking the equality of logarithmic derivatives in a 
chosen crossing point, whose location is immaterial for the result. The regular series solutions for 
the first and second equations of the even case are of the form yi{x) ~ Ai x''^ Si{x), i = 1, 2, where 
Ai are integration constants, Si{x) are power series in x with zeroth order term equal to unity and 
Pi are the positive indices 



1 1 

Pi = 



2 2(l-2e) 
1 1 



P2 = 



(1 - 2 ef {A~a^{4e + l))+4{l-2 e) j{j + 1)] 

1 1/2 
\2 „2 ' 



4(l + 2£)j(.7 + l)-(l + 2e)"a^ . (3.4) 



The asymptotic solutions for the two equations are of the form yi{x) — Bi exp[—Kix]x'^^Ti(x), 
1 = 1,2, where Bi are again integration constants, Ti{x) are power series in 1/x with zeroth order 
term equal to unity, while the two positive numbers Ki and the two indices Vi coincide for both 
equations and are equal to k and v, where 

1 



K = —\/~w (a^ w + 8) , 

^ _ w (8 + u;) (2 e + 1) + 16 
2 (4 + w) y^ur(8Tl?i(7) 

We finally remark that only the first terms of the series and asymptotic expansions have been 
determined analytically. The following ones have been calculated by numerical codes, and the 
number of terms to be taken has been chosen according to the following criterion: in the point where 
the initial conditions for the numerical integrations were assigned the result of the substitution of 
the solution into the differential equation divided by the solution itself had to be less that 10^^^. 
We also remark that tests have been made also for the arithmetic precision of the calculations and 
the number of meaningful figures has always been kept sufficiently high. 

Let us now discuss the solution of equation H3.2|l for the even ground state, j — 0, with e = 1/2, 
that presents a singularity at the positive point Xg = 2a/(4 + a'^w). We can study the series 
solutions in x^ and we realize that this singular point is in the case of the "limit cycle'^ of the Weyl 
classification [32] , independently of the value of w in the physical domain. This means that in the 
neighborhood of Xs the equation admits two finite solutions that can be used for matching the 
solution and its derivative coming from the origin to the solution and the corresponding derivative 
coming from infinity, forming therefore what in mathematics is referred to as a "classical solution" 
of the differential equation. Procedures of this type and also with more elaborate matching con- 
ditions have been occasionally considered, see e.g. |33| . but usually more as an investigation of 
mathematical possibilities, than under the real necessity of solving a specific problem. The first 
terms of the solution in the neighborhood of Xg are 

yg{x) =Aix- xg) (l + ^4^) +b(i+ ^-i---s)H\x-xAU ^^^^^ 



w + A J V w ^ A 
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The index in the origin is p = ^(1 + \/4^-3c?), while the parameters for the asymptotic solution 
are those given in (|3.5|) with e = 1/2. The results will be discussed in the next section. 

(b) The odd case 

From the remarks following H2.15|l . we see that the matrix for the odd system is obtained by 
substituting M ~ and fx = —2nie in H2.15|) . In the dimensionless variables 1)3. the explicit 
expression of the system we get is 

Tx " rh{x) + 2 h{x) = ° 

dx^^^""' xh{x)~2ae ^ x^^^""' ^ 2 x{xh{x)-2ae) 

I (h( ( A , 4j(j + l) 2 2 + 1) n 

:ry3W + ^ + - (4ae+ ^ -n-^) yiix) + -yz{x) — yi{x) = 

dx z V X 2ae — xn(x) / x xn[x) — 2ae 

d , , , 1 , 4j(j + l) , 4a^£2_^2;,2(^) 2 0-(j + l)y3(x) , 1 

+ 2 (^^7^ + x(xMx)-2a£)) ^) + ^ " 

(3.7) 

where —2-\- a^w/2 + a/a;. The system decouples only for = so that, for the moment, we 
will assume j > (remark that the triplet ground state has j = 1) . 

A superficial check of the coefficients shows the presence of an irrelevant singularity at positive 
values of x for e > 1/2. The situation, however, is not so simple and there actually exists a further 
hidden singularity that reveals itself when trying to integrate numerically the system in a direct 
way. This unexpected singularity becomes manifest when deducing the fourth order differential 
equation equivalent to the system (|3.7|l . Both the partial and the final results when obtaining this 
equation are very cumbersome: the final analytic expression have been calculated and managed 
only by means of a systematic use of computer algebra and cannot be reported here. We simply 
give the steps of the method we have followed to get the equation and that in mathematics is 
known as the prolongation method. First we isolate 2/2(2;) from the first equation and 2/4 (.t) from 
the third; we then substitute y2{x), its derivative and y^ix) in the second equation, obtaining a 
second order equation in yi{x) where yz{x) and its first derivative only appear. We then substitute 
1/4(0;), its derivative and 1/2(2^) in the fourth equation, obtaining a second order equation in 1/3(0;) 
where yi{x) and its first derivative only appear. We next consider the prolongations of this system: 
this means that we differentiate these two second order equations and substitute the second order 
derivatives obtained by the second order equation themselves. The third order equation for yi{x) 
, thus, contains 2/3 (x) and its first derivative only. We differentiate once more this equation and 
finally from this, the two third order and the two second order equations we can eliminate 2/3(2;) 
and its first, second and third derivatives, finally obtaining a fourth order differential equation for 

In the denominator of the coefficients of the fourth order equation we find a factor responsible 
for the appearance of the new singularity given by a root of the corresponding equation: 

-{a^w + 4) (16 + a^w {e'^ c? + 1) {c? w + 8)) x^ ~ 2 a(e^ u; (a^ u; + 8) (3 + 2 e) 
-{a^ w + Af{2e-i) + Z2e^ a^) x^ +4:a^ {-e'^ (4 e + 3) - 3 - 4 e (e^ - £ - l) 
+4e'^i(j + l))(a2 w + 4)x + 8a3 (4e''j(j + 1) - (2e + 1) (2e- 1)^ -e^ (2e+ i)^^) = 0. 

(3.8) 

For the typical values j = 1 and w = —0.5 it has a solution x > for e ~ 0.36014. Here also 
we stress that the presence of this singularity does not affect the perturbative approach to the 
hyperfine interaction, but must be considered if one is willing to integrate the odd equation as it 
stands. 

We have solved the fourth order equation for small enough values of e - so to prevent the 
presence of the additional singularity - as well as for e = 1/2. In both cases, taking into account 
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the due differences, the method is conceptually a variant of the double shooting procedure described 
in the previous paragraph. For small values of e we have to take care of the origin and infinity, 
that are the only two singular points. It turns out that in each one of them there exist two regular 
solutions that can provide the necessary initial conditions for starting the numerical integration. 
The spectral condition is now given by the matching of function, first, second and third derivatives 
in a fixed crossing point Xc, so to recover the "classical solution" of the equation. This amounts 
to the vanishing of the following determinant: 



det 



yoAXc) yoA^c) Voa,l{Xc) yoo,2{Xc) 



\ 



(i)^ ^ Wf \ (I) f \ (I) ^ ^ 
yo.i[Xc) Va^iK^c) Vj,i{xc) Voc,2\^c) 

(II)/ X (II)/ ^ (II)/ X (II)/ X 



Vyo.™^(^c) yo'2'ixc) y^'iixc) y^\2i^c) J 



,(III 



(III) 



,(III)^ 



= 



(3.9) 



where ?/o,ii yoo,i, (i — 1,2), are the two regular solutions coming from the origin and from the infinity 
respectively and the superscripts denote the order of the derivatives. As in the even case, the two 
acceptable series solutions in the neighborhood of the origin are of the form yi {x) = Ai x''^ St {x) , 
where, for j > and e < 1/2, 



P2 = 2 



2(1 - 2e) 



(1 - 2 e)2 (4 - a2 (4£ + 1)) - 4 (-1 + 2 e) j{j + 1) 



1/2 



-a^ (2e+l)2+4(2e+l)j(j + l) 



1/2 



(3.10) 



The search for the asymptotic solutions is a bit more delicate. Indeed they are still of the form 
yi{x) — Bi exp[— Kix] x^^Ti{x), i = 1, 2, and again ki = K2 = k, with 



-w (8 + w) 



1/2 



(3.11) 



while the indices are 



ly, = -1 



16 + a'^ w {8 + a'^ w) (2e + 1) 



2{A + a'^w) -w{8 + a'^w) 



1/2 



i^2 = I'l — 2 



(3.12) 



Although an integer difference of the indices could imply solutions of different type, the present 
case is the simplest one and a second non logarithmic solution is found. 

We now consider the odd ground state for e = 1/2. We have then to discuss the fourth order 
differential equation with j = 1 in the presence of a singularity located, almost independently of 
the value of w, around x ~ 0.0016478. The two regular asymptotic solutions are of the type already 
described and the corresponding parameters are obtained from (I3.11f) and H3.12|l with e = 1/2. 
The second solution is again non logarithmic. The behavior in the origin, however is here different. 
In fact we have one of the two regular solutions of the form yi{x) = x^'^~°'^ Si{x), but the second 
one, due to the fact that the singularity in the origin is irregular, must be searched in a more 
general form pi], and results in 



y2{x) = cxp 



4 ^/a{4^^ahu) 



(4 + a'^w) \fx 



x~^i^S2{^R) 



(3.13) 



where 6*2 (\/x) is a power series in y/x with zeroth order term equal to unity. Finally four regular 
solutions are found in a neighborhood of the singular point: their indices are 0,1,2,49/16, but 
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despite the integer differences all of them are non logarithmic. These solutions together with their 
first three derivatives are used to bridge the two solutions from the origin to the two solutions from 
infinity. From a numerical point, however, this is not very simple due to the critical sensitivity 
of the coefficients of the differential equation to tiny variations of the coordinate x. Indeed the 
ordinary integration codes present errors too large to be accepted, so that, in order not to lose 
in accuracy when integrating out of the singularity, we have chosen a mash of sufficiently close 
points (all of them obviously regular with respect to the differential equations): for each point we 
have constructed four regular series solutions with a number of terms sufficiently large to respect 
the accuracy requirements, and we have connected all these solutions by matching functions and 
derivatives up to order three. 

We finally discuss the spectral solution of the odd case with j = 0. We have already stated that 
now the system decouples and gives rise to a pair of second order differential equations. As in the 
even case, only one of these equations, namely the equation coming from the (2/1,2/3) subsystems, 
has a physical meaning. From H3.7|l we easily find 



X 



a 1 



yix) + - 1 + ^ — ' 5 TT^ TTT -ry(^) 



a^wx + 2a-\-8x Ax + a'^wx + 2a x aw + 2 



d 



dx 



- a'^ w {a"^ w + 8) a^{4e + l) a {a^ w + 4) {2 e + 1) 8ae 



16 4a;2 Ax + a'^ w) x + 2 a 



y{x) = 

(3.14) 



and we see that the coefficient of the first derivative has a singularity at the point x — ~2/{wa), 
independent of e. Since for the bound states we are studying w assumes negative values, the 
singularity is therefore located at a finite positive value of x and must be accounted for in the 
integration for any value of e. The situation is similar to what we have already seen and may 
briefly summarized as follows: the index for the acceptable series solution in the origin is p = 
— 1 + i ^4 — a^{4:e + 1). The two constants of the asymptotic solution are the same n as in H3.11|l 
and V = vi as given in (|3.12|l . In the neighborhood of the singular point the indices are 0,2 and 
there exist two non logarithmic solutions. 



4 Discussion of the results 

To begin the discussion of the results we report the values we have obtained in 27 for the levels of 
the pure Coulomb interacting system. The classification scheme we used in that paper was fit to 
describe the spectral terms for systems with variable mass ratio. For convenience in comparing our 
levels with the corresponding values existing in literature ^2 ^| and obtained by semi-classical 
expansions, we will adopt the commonly accepted classification of the Positronium levels, ^ . 

Our results for the relativistic two body equation with pure Coulomb interaction are as follows. 
For the ground states we have 



State 


W Coulomb 


State 


W Coulomb 


l^so 


-.4999950109 




-.4999950109 



For first excited states, the data are: 



^Although not necessary for the discussion, it is straighforward to recognize the following correspondence between 
El and IHj: (-h,0,/) ^ l^so, (-,!,/) ^ l^si, (+,0,11) 2^so, -> 2^pi, 2^pi, (-,0,1) 

23po, (-,!,//) ^ 23si, (-,2,7) ^ 23p2. 
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State 


W Coulomb 


State 


W Coulomb 




-.1249996884 


2^51 


-.1249996884 


2V1 


-.1250002427 


23pi 


-.1250005200 


23po 


-.1250007974 


23p2 


-.1249999654 



As explained in Appendix A, the levels accounting for the first order perturbative effects of the 
Breit terms are given by 



W Brcit = W Coulomb + 



1 dw{e) 

2 de 



e=0 



(4.1) 



The derivative of w in e = was calculated by computing the eigenvalues corresponding to 
e = 0.1 and e = 0.2 and then looking for the three point Lagrangian interpolation through them 
and through (e = 0, w = wcouiomb)- The results has been checked by repeating the same procedure 
for other values of e sufficiently close to the origin and the differences are at most of some unities 
on the last figure for the states s and even less for the states p. Even if we take the values we 
have calculated for e = 1/2 - and reported below - to construct the interpolation, we see that 
the results for the hyperfine ground levels differ only for some unities on the ninth figure. In the 
following table we summarize the data we have obtained. Remark that the calculations have been 
done with a number of figures sufficiently large to prevent the rounding errors. 



State 


We=0.1 


We=0.2 


dw/de\ e=o 


I'so 


-.5000008658 


-.5000024624 


-.7984077369- 10--* 


l^Sl 


-.4999955437 


-.4999953671 


-.8874896904-10-5 




-.1250005867 


-.1250009527 


-.1164413908-10-"^ 


2^si 


-.1249999215 


-.1250000659 


-.2774401032-10-5 


2V1 


-.1250003204 


-.1250002205 


-.1664091904-10-5 


2^pi 


-.1250007197 


-.1250008750 


-.2218611453-10-5 


2^Po 


-.1250012966 


-.1250017959 


-.4992271268-10-5 


2^P2 


-.1250000186 


-.1250000452 


-.6646165115-10-6 



The spectral values with e = 1/2 for I^sq and l^si are the following: 



State 


We=l/2 


State 


We=l/2 


l^So 


-.4999816915 


l^Sl 


-.4999905793 



The differential equations producing these last values present a singularity at finite values of the 
radial coordinate. This singularity is in fact non disturbing for the integration and one can "pass 
through" by imposing the matching conditions as explained in Section III. The results however 
confirms the Breit idea: indeed they present a qualitatively wrong configuration, with the singlet 
higher than the triplet. 

Let us now compare our results with known data. We remind the formula giving the first terms 
of the semi-classical expansion for the singlets 5 : 



w 



1 

'2^ 



2n4 V 16 



i + 



t)+o(«') 



(4.2) 
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It is well known that the l^si level has a value 

w^-- + —a'^ (4.3) 
2 96 ^ ' 

omitting the annihilation term, that contributes for an additional a^/2 to the singlet-triplet split- 
ting. From we take the approximations for the n = 2 triplets obtained from 10 . The results 
are 



and 



with 



w{2^p,) = w{2^po) + d{j) , u.(2Vi) = w{2^po) + 6' (4.5) 



S(l) ^ —o?^ ,5(2) = — a^ ,5' = —o?. (4.6) 

^ 160 ' ^ ^ 160 ' 24 ^ ' 

Finally w{2? s\) is deduced from the relation |5| 

w{2^s^) - u;(2iso) = \ (u'(l'si) - ii'(l'so)) (4.7) 
The comparison with our results is summarized in the following table: 



State 


^ num 


^ semi— classical 


I'so 


-.5000349313 


-i - =-.5000349462 


l^si 


-.4999994484 


-i + ^a2 =-.4999994453 


2^50 


-.1250055105 


-i - =-.1250055123 


2^31 


-.1250010756 


-i - T|ij.a2=-.1250010747 

8 1536 


2Vi 


-.1250010747 


-| - Y§ga2=-.1250010747 


2?px 


-.1250016293 


-| - yfga2^-.1250016294 




-.1250032935 


-i - liea^^ -.1250032935 


2?P2 


-.1250002977 


- 7||oa^=-.1250002982 



A couple of final comments on the results are in order. In the first place we observe that 
the pure Coulomb levels calculated in the semi-classical approximation or, equivalently, using the 
expansion in a differ from the levels calculated by the two body relativistic equation by a quantity 
of the same order of the approximation itself. In order to produce a concrete example we let e = 
in (|3.2|) and we use the "atomic variable" z — ^ax. Expanding in a up to the second order, we 
find 



dz^y^'^^y-^^^^^^2^ ' 4. 



^y{z)+[—^ + 2w+ + — ]y{z) (4.8) 



with regular solution 



y{z) = zV2 + (1/2) v/l-a%-(l/2) ^^w{a^w + S) z . 

1^1 i\ + ^ - , 1 + v/T^, v/-(«^- + 8).) . (4.9) 
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From the vanishing of the first argument of the hypergeometric we get the ground level 



w= ^(2 + 2^1- aA^2~2^l~a^~\^-\-^a^ + 0{a') (4.10) 

Comparing (|4.1U|) with our relativistic pure Coulomb result w — —.4999950109, that can be ap- 
proximated as w = — 1/2 + 3a^/32 + 0{a'^)^ we find a difference of about a^/A. On the contrary, 
when looking at the data presented in the final table, we see that the difference reduces to less than 
and becomes even more negligible for the levels with higher orbital angular momentum. 

A second observation concerns the degeneracy of the states 2'^si and 2^pi that is predicted from 
the perturbative expansion and that is completely confirmed from the numerical calculations. The 
existing difference of about 0.01 for the pure Coulomb interaction disappears when introduncing 
the Breit term. A similar phenomenon in the perturbative framework had been noticed in |14| 
with respect to the Lamb shift in the hydrogen atom, produced by a pure relativistic Coulomb 
interaction and reabsorbed by the presence of the magnetic term. 

To conclude, in this paper we have presented a relativistic calculation of the hyperfine structure 
of the Positronium, providing the theoretical and mathematical instruments to obtain the results. 
We have pursued our approach without any semi-classical approximation or expansion in the 
fine structure constant: the only due perturbative treatment has been reserved to the magnetic 
interaction term. Along this way, almost unexplored, we have proved some other aside facts. 
In particular we have discussed the nature and the properties of the singularities arising in the 
development and we have found that they bear no serious consequences neither in the integration 
of the wave equations, nor in their spectral behavior, but for lengthy technical complications: this, 
in a sense, can be considered an indirect test of the reliability of the approach to bound states 
through relativistic wave equations up to the quantum field theoretic corrections. We have also 
indicated possible applications of the method, that are now under investigation. 



A A perturbative expansion 

In this Appendix we prove the relationships between the derivative of the spectral values with 
respect to a parameter e and the perturbative expansion in that parameter. 
Consider the Hermitian operators H , Q and the sum 

K{e) = H + eQ (A.ll) 

Let U{£) be the unitary operator such that 

Kd{e) = U-\e)K{e)U{e) (A.12) 

is diagonal. Denoting by a dot the derivatives with respect to e and letting U = C^(0), U^^ = 
C/~^(0), U = C/(e)|e=Oi we expand equation (|A.12|I to the first order in e obtaining 

Kd + ekd = Hd + e([Ha, U-^U] + U-^QU^ (A.13) 

where Kj, — Kd{0), Kd = Kd{e)\^=Q. Obviously Kd = Hd = U^^HU and Kd are diagonal 
matrices. Moreover the diagonal elements of [Hd, U^^U] are vanishing, so that 

{Kdh^iU'^QU)u = {V,,QVi), [Hd,U-^U],j + {U'^QU),j^O, i^j, (A.14) 

where Vi is the i-th normalized eigenvector of H. Hence 

Kd{e)u = {Hdh + s{V,,QVi) + 0{e^) (A.15) 

and we recover the usual first order correction of the perturbative expansion. 

The procedure goes over to any order. Although in our framework we do not use anything but 
the first order, we want just to outline how the second order is also obtained. A straightforward 
computation gives 

Kd = [Hd, U-^U] + 2 U-^U[U-^U, Hd] + 2 [U-^QU, U-^U] (A.16) 
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where, in the usual notations, Kd = (P K d{e) / de'^]^^^ and the same for [7. Again the diagonal part 
of [Hd, U^^U] vanishes. Using the second equation of (jA.14|l and expanding the commutators, we 
find 

\{kd)u = Y.^U-^^),,{U-^QU),, + Y,iU-^QUh{U'^U)j^ - Y.^U-^^\,{U-^QU),, (A.17) 

3 3 

or equivalently, taking carefully into account the summation ranges, 

\{kd)u = ^(;7-iQC/)., (l/-i[/),, . (A.18) 



Observing that from the second equation of (|A.14p we have 



{U-'iJh ^ J^'^.^i". . (A.19) 
we find the usual second order contribution 

2 ^ ' ^ {Hdh - {Hd)n ^ ' 

It appears therefore the advantage achieved for a numerical evaluation of the perturbative correc- 
tions in our problem as previously explained. 
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